Topological and Dynamical Complexity of Random Neural Networks 
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Random neural networks are dynamical descriptions of randomly interconnected neural units. 
These show a phase transition to chaos as a disorder parameter is increased. The microscopic 
mechanisms underlying this phase transition are unknown, and similarly to spin-glasses, shall be 
fundamentally related to the behavior of the system. In this Letter we investigate the explosion 
of complexity arising near that phase transition. We show that the mean number of equilibria 
undergoes a sharp transition from one equilibrium to a very large number scaling exponentially 
with the dimension on the system. Near criticality, we compute the exponential rate of divergence, 
called topological complexity. Strikingly, we show that it behaves exactly as the maximal Lyapunov 
exponent, a classical measure of dynamical complexity. This relationship unravels a microscopic 
mechanism leading to chaos which we further demonstrate on a simpler class of disordered systems, 
suggesting a deep and underexplored link between topological and dynamical complexity. 

PACS numbers: 87.18.Tt, 05.10.-a 87.19.11, 87.18.Sn, 87.18.Nq 



Heterogeneity of interconnections is a crucial property 
to understand the behavior of realistic networks that 
arise in the modeling of physical, biological or social 
complex systems. Among these, a paramount example 
is given by neuronal networks of the brain. In these sys- 
tems, synaptic connections display characteristic disorder 
properties [l], Q resulting from development and learn- 
ing. Taking into account this heterogeneity seems now 
essential, as experimental studies of neuronal tissues have 
shown that the degree of disorder in the connections sig- 
nificantly impacts the input-output function, rhythmicity 
and synchrony, effects that can be related to transitions 
between physiological and pathological behaviors. As an 
example, Aradi and Soltesz [3] have shown that rats sub- 
ject to febrile seizure present the same average synaptic 
properties but increased variance. 

These properties are reminiscent of disordered phys- 
ical systems such as spin-glasses. The relationship be- 
tween disorder and qualitative behaviors has been thor- 
oughly studied within theoretical frameworks such as the 
Sherrington-Kirkpatrick model 1-1] describing the behav- 
ior of binary variables interacting through a random con- 
nectivity matrix. In these models, estimating the number 
of metastable states, deemed to be deeply related to the 
transition to chaos, remains an important endeavor (EJ. 
In the context of nervous system modeling, random neu- 
ral networks 0,0 constitute a prominent class of models, 
which describe the evolution of the activity of a neuron 
i in a n-neurons network through the randomly coupled 
system of ordinary differential equations: 

n 

X{ = — X\ + ^ ^ JijS(Xj) (1) 
j=l 

where Jij are independent centered Gaussian random 
variables of variance a 2 /n representing the synaptic con- 
nectivity coefficients between neuron i and j, and S is 



an odd sigmoid function with maximal unit slope at the 
origin (S'(x) < S'(0) = 1) representing the synaptic non- 
linearity. 

The behavior of system (pQ) has been analyzed in the 
asymptotic regime of infinite population size n — >• oo Q 
and displays a phase transition for a critical value of the 
disorder at a = 1: for a < 1 all the trajectories are at- 
tracted to the trivial equilibrium x = and for a > 1 the 
trajectories have chaotic dynamics. This generic phase 
transition has been numerically observed in a number 
of situations in more realistic models involving multiple 
populations and excitable dynamics [8j . That phase tran- 
sition, and the chaotic regime beyond the edge of chaos, 
appear particularly relevant to understand the computa- 
tional capabilities of neuronal networks. In particular, 
information processing capacity was characterized as op- 
timal at the edge of criticality 0, [To| , and such random 
neural networks are in particular used in recent machine 
learning algorithms [ll|. Moreover, the question of crit- 
icality has been widely debated in the theoretical neu- 
roscience community and beyond [12l-ll4j and it seems 
that through a number of different mechanisms among 
which plasticity of synapses, networks may tend to be 
naturally poised near criticality. In the random neural 
network, the microscopic mechanisms underpinning this 
phase transition remains unclear. It is hence of great in- 
terest to dissect precisely the behavior of such systems 
at the edge of chaos. 

This question of the topological modifications of the 
phase space arising at the edge of chaos is precisely the 
question we shall address in this Letter. More precisely, 
we will characterize the averaged number of equilibria in 
the random neural network (pQ). To this end, we will de- 
velop upon the theories of random matrices [15J , random 
geometry and Gaussian fields [16|- 

We denote by A n (a) the random number of equilibria 
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(depending on the realization of the matrix J = (Jij)). 
These are the solutions of the system: 



JijS(xj) 



(2) 



For a < 1, consistently with the mean-field analysis, 
we first show that E[A n (cr)] — > 1 when n — >> oo. The 
proof proceeds by showing that the system is contracting 
on the whole space R n . To this purpose, one needs to 
characterize the eigenvalues of the random Jacobian ma- 
trix — I + J.A(S"(x)) where I is the identity matrix, and 
A(5"(x)) is the diagonal matrix with elements S'(xi). 
The matrix J.A(S"(x)) is a centered Gaussian matrix 
with independent components, and each column has a 
distinct standard deviation given by cr 2 S' (xi) 2 /n. Raj an 
and Abbott [ItJ provide the system of equations satis- 
fied by the squared modulus of the eigenvalues of such 
random matrices, and solve these when considering only 
two different variances. Their methodology readily gen- 
eralizes to our problem, and elementary algebraic manip- 
ulations (see 3, [HI) show that the spectral density has 
support in the disc of radius — Yl7=i S'( x i) 2 f° r large 
n. In our case, S'(x) < 1, and hence all eigenvalues of 
the matrix — I + J • A(5"(x)) have a negative real part in 
the limit n — >> oo. This implies global contraction of the 
dynamics, ensuring the fact that the trivial equilibrium 
x = is the unique equilibrium for a < 1 (Banach fixed 
point theorem) and its global stability. 

For a > 1 the situation is more complex. Similarly 
to the phase transition in spin glasses, the behavior of 
E[j4 n (cr)] is likely to scale as exp(nC(a)) where C(a) 
is the topological complexity. Our starting point is to 
observe that the fixed point equations (|2j) can be viewed 
as zero crossings of a Gaussian field indexed by x G M n , 
and therefore we may apply the Kac-Rice [16j formula in 
order to estimate the number of solutions: 



E[A n (<r) 



dxE 



det(-I + J.A(S"(x)))| 



X(J (-x + J.5(x)) 



Recent studies [2CH22] have used the Kac-Rice formula 
to estimate the mean number of critical points in ran- 
dom energy landscapes arising from Hamiltonian spin- 
glass with strong symmetry properties. In our case, the 
situation is substantially different: (i) there is no under- 
lying energy landscape since the system (pQ) is not Hamil- 
tonian and (ii) symmetry properties are relatively weak 
and do not enable the same kind of reduction method 
developed in [2CH22|. A major technical difficulty is the 
fact that, because the lack of symmetry in our system, 
one needs to deal with determinant of random matrices 
with columns of non-identical variances, for which the 
spectral density is unknown. 

However, near criticality, for a = 1 + s with 0<£<1, 
we can obtain a first order estimate of the number of 



equilibria taking advantage of the fact that all equilibria 
remain close to x = 0. More precisely, we first show 
that with an arbitrarily high probability 1 — £(n, e) with 
£(n,e) — » as n —> oo, all equilibria belong to a ball 
B p u\ centered at x = of radius p(e) which tends to as 
e — >• 0. This is a consequence of the spectral analysis of 
the random matrix — I + J • A(5"(x)), whose eigenvalues 
all have negative real parts outside of B p ( e ) for large n. 
The property that fixed points remain in a small ball 
around zero is non-trivial. The proof proceeds by defining 
p(e) the unique positive solution of the scalar equation 
x/a = S(x), which is clearly arbitrarily small when e — >• 
0, and the smooth modified sigmoid function S £jTJ 



x/a 
S(x) 



\x\ < p(e) 
\x\ > p(e)+r] 



(the small interval [p(e),p(e) + rj] allows to define a 
smooth continuation). Because of the properties of the 
sigmoidal function, and in particular the fact that the 
differential is decreasing for x > (and increasing for 
x < 0), the same argument as used in the case a < 1 
applied to system (pQ) defined with the sigmoid S £iTJ 
(termed modified system) ensures that all eigenvalues of 
— I + J • A(S' £ ^(x)) have negative real part and hence 
that the unique fixed point of the modified system is 0. 
In particular, there is no fixed point outside the ball of 
radius p(e). In that region the original and modified sys- 
tem are identical, implying that the only possible fixed 
points of the original system are contained in the ball of 
radius p(e). 

Therefore, one can split the expectation according to 
whether |x| < p(e) or not, yielding: 



E[A 



dxE 



det(-I + J.A(S"(x)))| 



X(5 (-x + J.S(x)) +0(f(n,e)) 



Moreover, thanks to the differentiability of the determi- 
nant operator, we know that within the ball B p ( e ), the 
integrand is equal to | det(— 1 + J)| +0(p(e)), eventually 
yielding: 

E [A n (a)} = E [| det(-I + J)|] + 0(p(e) + e)) 

To evaluate this formula, we first compute the logarithm 
of the determinant: 

-log|det(-I + J)| = - lo sl A - !| 

\esp{j) 

where sp(J) denotes the spectrum of J, which in the large 
n limit is uniformly distributed in the disc of radius a 
(23| . Using this property one obtains in the large n limit 



1 



log I det(-I + J) I = c(cr) + R(n) 
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and c(a) := / c log \z — l\fi a (dz) with fi a (dz) the uniform 
distribution on the disk of radius a = 1 + e and R(n) is 
the finite-size error associated to the convergence to the 
circular law. We have: 



■ logE [| det(-I + J)|] = c{&) + - logE 



iR(n) 



■ (3) 



This ensures that ^ logE [A n (cr)] is arbitrarily close c(a) 
when e — >> oo and n oo. 

We are hence left computing c(a). Since log \z\ is har- 
monic, we can show that for a, 6 > : 



\og\a-be w \dO = 27rlog(min(a, b)) 



ensuring that c(a) = for a < 1, which is consistent with 
our previous analysis of this case, and 



c(a) = 



1 



2nr log(r) 



log(a) + 



1 



1 



(4) 



for a > 1. For a close to 1 + , we conclude c(a) ~ (a — l) 2 . 
Therefore we have shown that: 



E[A n (<7) 



n(cr-l) 2 



up to multiplicative polynomial factors. 

This combinatorial explosion of the number of equi- 
libria is the hallmark of the accumulation of bifurcations 
[29| in a the neighborhood of the critical parameter value 
(7 = 1. Moreover, this result provides a possible topolog- 
ical explanation for the emergence of chaos. For a > 1, 
the phase space is heavily mined with equilibrium points, 
most of which are saddles due to the spectral properties 
of the Jacobian matrices as discussed above. Typical 
trajectories evolving in this landscape will hence wander 
from the vicinity of one point to another following sta- 
ble manifolds. The presence of a large number of saddle 
points induces a very high sensitivity to perturbations, 
distinctive feature of chaos. 

The classical characterization of chaos relies on the 
evaluation of the maximal Lyapunov exponent of the 
trajectories quantifying the dynamical complexity. This 
quantity is defined as follows. Applying an infinitesi- 
mal external perturbation X{ — >• X{ + 5x® on neuron i 
at time to induces a change on all neurons at subse- 
quent times Xj (t) — )• Xj (t) + Sxji (t) defining a suscepti- 
bility matrix ^^(to + r, to) = Sxji(to + r)/5x^. From 
this matrix we shall define an averaged susceptibility 



\I/ 2 (t) lim 



t - 



£i, ? -*ii(*0+T,t ) 



ally define the Lyapunov exponent by: 



and eventu- 



lim 



log(* 2 (£)) 



This quantity was analyzed in [6] using a spectral de- 
composition based on the one-dimensional Schrodinger 



equation. Near criticality, the decomposition dramati- 
cally simplifies and yields for a ~ 1+: 



A ~ (cr 



i) 2 . 



We hence conclude that the topological and dynamical 
complexities have the same behavior at the edge of chaos. 
More surprising is that using the spectral decomposition 
of the Schrodinger equation in the limit a ^> 1, Sompolin- 
sky and collaborators show that A diverges as log(cr), 
precisely as the complexity of the system given in for- 
mula (|4]), although our analysis rigorously only applies 
for a close to 1. This further strengthens our microscopic 
interpretation of the emergence of chaos in relationship 
with the number of saddles. 

If this interpretation is the actual phenomenon aris- 
ing in random neural networks, then simpler, lower- 
dimensional dynamical systems shall display a similar re- 
lationship between the number of unstable equilibria and 
the Lyapunov exponent. Probably the simplest model 
one could think of is what we call the fakir bed dynam- 
ics, corresponding to the movement of a particle in a 
two-dimensional complex landscape with k unstable fixed 
points. More precisely, we now consider a particle con- 
fined in a compact subset of R 2 , with close to the origin 
a fixed number k of Gaussian hills (corresponding to the 
presence of k unstable fixed points) randomly chosen in 
space. Trajectories of particles soon get chaotic as k in- 
creases (see Fig. [2(b)), and we numerically compute the 
maximal Lyapunov exponent of the trajectories. Since 
the landscape is probabilistic due to the choice of the 
location of hills, we compute the averaged maximal Lya- 
punov exponent across 100 independent realizations of 
the process, and plot it against the logarithm of the num- 
ber of equilibria. The corresponding curve (Fig. QJa)), 
indeed displays an increasing profile well approximated 
by a linear curve: a clear relationship again emerges be- 
tween number of unstable fixed points and Lyapunov ex- 
ponents, supporting our interpretation related to the ran- 
dom neural network. 





(a) Lyapunov exponent vs 
number of fixed point 



(b) Fakir bed trakectories 



FIG. 1. The Fakir bed dynamics, (a): blue: Lyapunov expo- 
nent averaged over 100 realizations of the Fakir's dynamics vs 
number of fixed points (in semilogarithmic axes) and linear 
regression, (b): Red and Green: two trajectories with very 
close initial conditions. 

In this letter, motivated by the analysis of fine micro- 
scopic phenomena arising at the edge of chaos in random 



neural networks, we have shown that the average num- 
ber of equilibria scales exponentially with the system size, 
with an exponential coefficient proportional to the Lya- 
punov exponent. This property, relying in part on the 
spectral theory of random matrices, readily inherits uni- 
versality properties of the circular law [25| , and therefore 
is valid for a large class of independent couplings beyond 
Gaussian. 

The complexity was shown to vary smoothly with the 
disorder parameter with a critical exponent 2, larger 
than 1, ensuring a form of structural robustness in the 
neighborhood of the phase transition. This study is the 
first application to random neural networks of recent 
methods used for counting the number of metastable 
equilibria in spin glasses. From the theoretical view- 
point, we extended that approach to out of equilibrium, 
non-Hamiltonian systems at zero temperature (singu- 
lar points of a vector field). Incidentally, the relation- 
ship found between topological and dynamical complex- 
ity highlighted what we conjecture to be a deep corre- 
spondence. Numerous open questions and perspectives 
have emerged from this study, among which the estima- 
tion of the distribution of equilibria and their number 
beyond the edge of criticality, requiring significant ad- 
vances in the analysis of random matrix determinants, 
or pursuing the exploration, in line with [26l-[28|. of the 
relationship between topological and dynamical complex- 
ity with other measures such as the fractal dimensions of 
chaotic attractors. 
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